# for Vlasi look for ##  comments forthe family ##
# last change 2-4-11
# garmaFit() is a function based on the Garch11fit() of the fGarch pacakge 
# it uses a non linear maximation routine to fit a GARMA(1,1) model
# it is a usefull application of the recursive use of filter() 
# TO DO
# 
# 1) Can generalised for more general ARMA model arma(p,q)? yes  OK
# 2) for more general distributions yes OK
# 3) needs se error for te initial parametets so we can get lower and upper bounds: maybe we dont need this set to -Inf nad Inf
# 4) residuals() OK
# 5) test different link functions for g2(y*) 
# 6) fix to allows 3 and 4 parameter distributions
# 7) fitted()   OK
# 8) deviance() OK
# 9) AIC()      OK
#10) vcov()     OK
#11) summary()  OK
#12) plot()     OK
#13) print()    ok
#14) coef()     ok
#15) wp()
# 16) fixed parameters
source("/Users/stasinom/Documents/gamlss/projects/TIME-SERIES/createLags.R")
#---------------------------------------------------------------------------------------
#---------------------------------------------------------------------------------------
require(gamlss)
#---------------------------------------------------------------------------------------
#---------------------------------------------------------------------------------------
garmaFit <- function(formula = formula(data), 
		               order = c(0,0),
		             weights = NULL,
					    data = sys.parent(), 
					  family = NO(),  ##  declaring the family ##
					   alpha = 0.1,
				   phi.start = NULL,
				 theta.start = NULL, 
					    tail = max(order),
					 control = list())
{
#---------------------------------------------------
# local functions
#--------------------------------------------------------------------------------------
# this is to replicate rqres within gamlss enviroment DS Friday, March 31, 2006 at 10:30
# it is used as in gamlss()
	rqres <- function (pfun = "pNO", 
			type = c("Continuous", "Discrete", "Mixed"),
			censored = NULL,  
			ymin = NULL, 
			mass.p = NULL, 
			prob.mp = NULL,
			y = y,
			... )
	{ }
	body(rqres) <-  eval(quote(body(rqres)), envir = getNamespace("gamlss"))
	##---------------------------------------------------------------------------------------
	##---------------------------------------------------------------------------------------
# the hessian	
#---------
	hessian<-function(par)
	{
		npar <- length(par)
		epsilon <- 0.0001 * par
		Hessian = matrix(0, ncol = npar, nrow = npar)
		for (i in 1:npar) 
		{
			for (j in 1:npar) 
			{
				x1. <- x2. <- x3. <- x4. <- par
				x1.[i] <- x1.[i] + epsilon[i]; x1.[j] <- x1.[j] + epsilon[j] 
				x2.[i] <- x2.[i] + epsilon[i]; x2.[j] <- x2.[j] - epsilon[j]
				x3.[i] <- x3.[i] - epsilon[i]; x3.[j] <- x3.[j] + epsilon[j]
				x4.[i] <- x4.[i] - epsilon[i]; x4.[j] <- x4.[j] - epsilon[j]
				Hessian[i, j] <- (garmaLL(x1.)-garmaLL(x2.)-garmaLL(x3.)+garmaLL(x4.))/(4*epsilon[i]*epsilon[j])
			}
		}
		Hessian
	}
#----------
	garmaLL <- function(parm,  save=FALSE) 
	{	
		           beta <- parm[1:l.beta]	
		             Xb <- X%*%beta
		switch(case, 
				{ # case 1
			        start  <- l.beta + 1
				    finish <- start + (order[1]-1)
				       phi <-  parm[start:finish]
				  finish   <- finish + 1
				    g2ylag <- createLags(g2y, order[1])
				     Xblag <- createLags(Xb, order[1])
				    g2y_Xb <- g2ylag-Xblag
				       eta <- Xb +  g2y_Xb%*%phi
				        mu <- ifelse(is.na(eta), alpha, fam$mu.linkinv(eta)) 
				},
				{ # case 2
					
					start  <- l.beta + 1
					finish <- start+ (order[2]-1)
					 theta <- parm[start:finish]
				    finish <- finish + 1
					g2ylag <- createLags(g2y, order[2])     
					  some <- Xb+g2ylag%*%theta
					   eta <- filter(some, -theta, "r", init =rep(0,order[2]))
	#				   eta <- try(filter(some, -theta, "r", init =rep(0,order[2])))
	#				   if (any(class(eta)%in%"try-error"))
    #             { 
    #               browser()
    #
    #             } 
					    mu <- ifelse(is.na(eta), alpha, fam$mu.linkinv(eta)) 	
				},
				{ # case 3
					#print(parm)
					#ƒif(any(is.na(parm))) browser()
					  start <- l.beta + 1
					finish1 <- start + (order[1]-1)
					finish2 <- finish1+1+ (order[2]-1)
					   phi  <- parm[start:finish1]
					  theta <- parm[(finish1+1):finish2]
					 finish <- finish2+1    
				   g2ylagar <- createLags(g2y, order[1])    
				   g2ylagma <- createLags(g2y, order[2])
				    Xblagar <- createLags(Xb, order[1])
				    Xblagma <- createLags(Xb, order[2])
				  	 g2y_Xb <-  g2ylagar-Xblagar
				  	#print(phi)
				  	# print(theta)
				  	 #print(head(cbind( g2y_Xb, g2ylagar, g2ylagma, Xblagar, Xblagma)))
					   some <- Xb+g2y_Xb%*%phi+g2ylagma%*%theta
					  # print(head(some))					  
					 eta <- filter(some, -theta, "r", init =rep(0, order[2])) 
			#tryCatch(eta <- filter(some, -theta, "r", init =rep(0, order[2])), finally=print(head(eta)))
					   mu <- ifelse(is.na(eta), alpha, fam$mu.linkinv(eta)) 		  
				}
		       )
			   switch(nopar, ##  this is where swiches between different family with different distributions ##
					   {# one parameter 
						   llh <- -sum(w*PDF(y, mu=mu,  log=TRUE), na.rm = TRUE)
					   },
					   {# two  parameters 
						 sigma <- parm[finish] 
						   llh <-  -sum(w*PDF(y, mu=mu, sigma=sigma, log=TRUE), na.rm = TRUE)
					   },
					   {# three parameters
						 sigma <- parm[finish] 
						    nu <- parm[finish+1]
						   llh <- -sum(w*PDF(y, mu=mu, sigma=sigma, nu=nu, log=TRUE), na.rm = TRUE)
						   
					   },
					   {# four parameters
						 sigma <- parm[finish] 
						    nu <- parm[finish+1]
						   tau <-parm[finish+2]
						   llh <-  -sum(w*PDF(y, mu=mu, sigma=sigma, nu=nu, tau=tau, log=TRUE),  na.rm = TRUE)
					   }
			         )
		if (save) return(list(lik=llh, mu=mu))
	    #if(is.na(llh)) browser()
		llh 
	}	
#---------------------------------------------------------
#---------------------------------------------------------
# the main function starts here
	garmacall <- match.call()  #   the function call
## which case of garma?
##  order=c(0,0)     case 0
##  order=c(>=,0)    case 1
##  order=c(0,>0)    case 2
##  order=c(>0,>0)   case 3
      if (order[1]<=0&&order[2]<=0)  case <- 0 
      if (order[1]> 0&&order[2]<=0)  case <- 1 
	  if (order[1]<=0&&order[2] >0)  case <- 2	
	  if (order[1]> 0&&order[2]>0 )  case <- 3
## starting values for beta
## fitting a gamlss model
## possibly with weights ?????
     m0 <- gamlss(formula, family=family, data=data, trace=FALSE) #with(data, gamlss(y ~ x2 + x3 + x4 + x5, family=family,  trace=FALSE)
     cat("deviance of linear model= ", deviance(m0),"\n")
     if (case==0) return(m0) # stop if case 0
# get the initial betas
       beta <- coef(m0)
	   # here we need the se so we can create lower and upper bounds???
	#se.coef <-  c(0.09662, 0.13091, 0.13968, 0.13302, 0.13437) 
	#   se.coef <- vcov(m0, "se") NEEDS ATTENTION
	 l.beta <- length(beta)
	#beta.se <- [1:l.beta]
          y <- m0$y # 
	      X <- m0$mu.x # the X matrix
	    fam <- as.gamlss.family(family)
	#if (!is.null(data)) {attach(data); on.exit(detach(data))}
   ##  the next 6 lines are needed for the  family ##
      fname <- fam$family[[1]] 
       dfun <- paste("d",fname,sep="")
       pfun <- paste("p",fname,sep="")
	    PDF <- eval(parse(text=dfun))
	    CDF <- eval(parse(text=pfun))
      nopar <- fam$nopar
     	  N <- length(y)
	# depending on the link function
	# I am not sure I cover all possibilities what about own???
	# also the logit needs testing ?????
	     if (fam$mu.link=="identity") ystar <- y 
	     if (fam$mu.link=="log")      ystar <- pmax(y,alpha)
	     if (fam$mu.link=="logit")    ystar <- ifelse(y==0, alpha, ifelse( y==1, 1-alpha, y) )
	     if (fam$mu.link=="inverse")  ystar <- y
	    g2y <- fam$mu.linkfun(ystar) 
	# I will need here ar and ma 
    tailoff <- max(order)
	  mtail <- max(tailoff, tail)
	      w <- if(is.null(weights))    rep(1, N) else weights
	          if(any(w < 0)) stop("negative weights not allowed") # 
 w[1:mtail] <- 0	# weight out the tail  
# get the initial values if are not set by the user
##  you will need  something like this for the general family probable you need for nu and tau only##
	if ("mu"%in%names(fam$parameters))
	{
		params <- c(beta = beta) 
   lowerBounds <- c(beta = rep(-Inf, l.beta)) # those bounds need checking ???
   upperBounds <- c(beta = rep( Inf, l.beta))
		 switch(case, 
				 {
				 	
				 	if (length(phi.start)!=order[1]&&!is.null(phi.start)) 
				 	        stop("phi.start should have ", order[1], " elements \n" )
				      phi <- if (is.null(phi.start)) runif(order[1]) else phi.start           
			       params <- c(params,       phi = phi)
			  lowerBounds <- c(lowerBounds,  phi = rep(-1,order[1]))
			  upperBounds <- c(upperBounds,  phi = rep(1, order[1])) 
			      },
		          {
		          	if (length(theta.start)!=order[2]&&!is.null(theta.start)) 
				 	        stop("theta.start should have ", order[2], " elements \n" )
				    theta <- if (is.null(theta.start)) runif(order[2]) else theta.start  
			       params <- c(params,      theta= theta)
			  lowerBounds <- c(lowerBounds, theta = rep(-1,order[2]))
		      upperBounds <- c(upperBounds, theta = rep(1, order[2])) 	  
			      },
				  {
				  	if (length(phi.start)!=order[1]&&!is.null(phi.start)) 
				 	        stop("phi.start should have ", order[1], " elements \n" )
				 	if (length(theta.start)!=order[2]&&!is.null(theta.start)) 
				 	        stop("theta.start should have ", order[2], " elements \n" )
				 	   phi <- if (is.null(phi.start)) runif(order[1]) else phi.start 
				 	 theta <- if (is.null(theta.start)) runif(order[2]) else theta.start           
				   params <- c(params,       phi = phi,   theta = theta)
		      lowerBounds <- c(lowerBounds,  phi = rep(-1,order[1]),  theta = rep(-1,order[2]))
		      upperBounds <- c(upperBounds,  phi = rep(1, order[1]),  theta = rep(1,order[2])) 	  
				  }
		       )
    }
	if ("sigma"%in%names(fam$parameters))
	{
		     sigma <- fitted(m0, "sigma")[1] 
	  names(sigma) <- ""
		#sigma.coef <- coef(m0, "sigma")
		#  sigma.se <- se.coef[l.beta+1]
		# This needs fixxing  
		#lowersigma <-  fam$smgma.linkinv()
		params <- c(params, sigma = sigma )
		lowerBounds <- c(lowerBounds, sigma =  0.001)
		upperBounds <- c(upperBounds, sigma =  1e10) 
	}
	if ("nu"%in%names(fam$parameters))  ##  probably ypu need something similar ##
	{
		nu <- fitted(m0, "nu")[1] 
		names(nu) <- ""
		#sigma.coef <- coef(m0, "sigma")
		#  sigma.se <- se.coef[l.beta+1]
		# This needs fixxing  
		#lowersigma <-  fam$smgma.linkinv()
		params <- c(params, nu = nu )
		lowerBounds <- c(lowerBounds, nu = -Inf)# I have to think about this ?
		upperBounds <- c(upperBounds, nu =  Inf) 
	}
	if ("tau"%in%names(fam$parameters)) ##  probably ypu need something similar ##
	{  
		tau <- fitted(m0, "tau")[1] 
		names(tau) <- ""
		params <- c(params, tau = tau )
		lowerBounds <- c(lowerBounds, tau = -Inf)# I have to think about this ?
		upperBounds <- c(upperBounds, tau =  Inf)
	} 
	     # S <- 1e-6
	
#  Estimate Parameters 
##  now we have to make sure that the likelihood can cope with different famillies so see garmaLL ##
#	fit <- nlminb(start = params, objective = garmaLL, lower = lowerBounds,  upper = upperBounds, control = control)
 				   fit <- try( nlminb(start = params, objective = garmaLL, lower = lowerBounds,  upper = upperBounds, control = control))
					   if (any(class(fit)%in%"try-error"))
                 { 
                 cat("OOPS IT FAILED: let us try once more \n")	
                  fit <- do.call("garmaFit", args=as.list(garmacall[-1]))
                  return(fit)
                 } 
   
    
    # fit <-  optim(params, garmaLL, method = c("L-BFGS-B"),lower = lowerBounds, upper = upperBounds, control = list(), hessian = FALSE)
   fit1 <- garmaLL(fit$par, save=TRUE)  
 df.fit <- length(fit$par)
Hessian <- hessian(fit$par)   
# Step 6: Create and Print Summary Report:
#          se.coef <- sqrt(diag(solve(Hessian)))
#             tval <- fit$par/se.coef
#          matcoef <- cbind(fit$par, se.coef, tval, 2*(1-pnorm(abs(tval))))
#dimnames(matcoef) <- list(names(tval), c(" Estimate",
#" Std. Error", " t value", "Pr(>|t|)"))
#cat("\nCoefficient(s):\n")
#printCoefmat(matcoef, digits = 6, signif.stars = TRUE)
# preparing for output
cat("deviance of  garma model= ", 2*fit$objective,"\n")
if (2*fit$objective>deviance(m0)) warning("There is a problrm with the fitted GARMA model here \n")
out <-list(family = fam$family,  
       parameters = as.character(names(fam$par)), 
		     type = fam$type, 
		     call = garmacall, 
		        y = y,
	      weights = w,
       G.deviance = 2*fit$objective, 
	       	    N = N,  
	       df.fit = length(fit$par), 
      df.residual = N-df.fit, 
	 	      aic = 2*fit$objective+2*df.fit,
		      sbc = 2*fit$objective+log(N)*df.fit,
	       method = "nlminb",
		     vcov = solve(Hessian),
		     coef = fit$par) 
if (   "mu"%in%names(fam$parameters))      	
{		
     		      mu <- as.vector(fit1$mu)
	    # mu[1:mtail] <- NA
 		      out$mu <- mu
			  if (fam$nopar>1) 
		       {
				   final <-which(names(fit$par)=="sigma")-1
     out$mu.coefficients <- fit$par[1:final]
		       }
		      else 
			  {  
				  out$mu.coefficients <- fit$par
			  }
	  }
if ("sigma"%in%names(fam$parameters)) 
      {
	  	        sigma <- fit$par[ "sigma"]
	        out$sigma <- sigma
out$sigma.coefficients<- fam$sigma.linkinv(sigma)
	  }
if (   "nu"%in%names(fam$parameters))    
	  {
		         nu <- fit$par[ "nu"]
	         out$nu <- nu
out$nu.coefficients <- fam$nu.linkinv(nu)		 
	  }
if (  "tau"%in%names(fam$parameters))
      {
	             tau <- fit$par[ "tau"]
	         out$tau <- tau
out$tau.coefficients <- fam$tau.linkinv(tau)	
      }
	out$residuals <- eval(fam$rqres)
	    out$rqres <- fam$rqres
	class(out) <-c("garma", "gamlss")
	out
}
######################################################################################
#-------------------------------------------------------------------------------------
# garmafit finish here
#-------------------------------------------------------------------------------------
######################################################################################
vcov.garma<-function (object, type = c("vcov", "cor", "se",  "all"), ... ) 
{
	type <- match.arg(type)
	switch(type,
			"se"={x <- sqrt(diag(object$vcov))},
			"cor"={x <- cov2cor(object$vcov)},
			"vcov"={x<- object$vcov},
			"all"={x <-list(vcov=object$vcov, cor=cov2cor(object$vcov), se=sqrt(diag(object$vcov)))}
	)
	x
}
######################################################################################
######################################################################################
fitted.garma<-function (object, what = c("mu", "sigma", "nu", "tau"), ... ) 
{
	what <- match.arg(what)
	if (! what%in%object$par) stop(paste(what,"is not a parameter in the gamlss object","\n"))
	x <- if (what=="mu") as.vector(object$mu)
	     else rep(object[[what]], object$N)
	x
}
######################################################################################
######################################################################################
print.garma <- function (x, digits = max(3, getOption("digits") - 3), ...) 
{
    cat("\nFamily: ", deparse(x$family), "\n")
    cat("Fitting method:", deparse(x$method), "\n")
    cat("\nCall: ", deparse(x$call), "\n", fill = TRUE)
    cat("Mu Coefficients")
    if (is.character(co <- x$contrasts)) 
        cat("  [contrasts: ", apply(cbind(names(co), co), 1, 
            paste, collapse = "="), "]")
    cat(":\n")
    if ("mu" %in% x$parameters) {
        print.default(format(coef(x, "mu"), digits = digits), 
            print.gap = 2, quote = FALSE)
    }
    if ("sigma" %in% x$parameters) {
        cat("Sigma Coefficients:\n")
        print.default(format(coef(x, "sigma"), digits = digits), 
            print.gap = 2, quote = FALSE)
    }
    if ("nu" %in% x$parameters) {
        cat("Nu Coefficients:\n")
        print.default(format(coef(x, "nu"), digits = digits), 
            print.gap = 2, quote = FALSE)
    }
    if ("tau" %in% x$parameters) {
        cat("Tau Coefficients:\n")
        print.default(format(coef(x, "tau"), digits = digits), 
            print.gap = 2, quote = FALSE)
    }
    cat("\n Degrees of Freedom for the fit:", x$df.fit, "Residual Deg. of Freedom  ", 
        x$df.residual, "\n")
    cat("Global Deviance:    ", format(signif(x$G.deviance)), 
        "\n            AIC:    ", format(signif(x$aic)), "\n            SBC:    ", 
        format(signif(x$sbc)), "\n")
    invisible(x)
}
######################################################################################
######################################################################################
coef.garma <- function (object, what = c("mu", "sigma", "nu", "tau"), ...) 
{
    what <- match.arg(what)
    if (!what %in% object$par) 
        stop(paste(what, "is not a parameter in the object", 
            "\n"))
    x <- object[[paste(what, "coefficients", sep = ".")]]
    x
}#
#####################################################################################
#####################################################################################
summary.garma <- function (object, digits = max(3, getOption("digits") - 3), ...) 
{
    cat("\nFamily: ", deparse(object$family), "\n")
    cat("Fitting method:", deparse(object$method), "\n")
    cat("\nCall: ", deparse(object$call), "\n", fill = TRUE)
   # cat("Mu Coefficients")
   # if (is.character(co <- object$contrasts)) 
   # cat("  [contrasts: ", apply(cbind(names(co), co), 1, 
   #         paste, collapse = "="), "]")  
   # Step 6: Create and Print Summary Report:
          se.coef <- sqrt(diag(object$vcov))
             tval <- object$coef/se.coef
          matcoef <- cbind(object$coef, se.coef, tval, 2*(1-pnorm(abs(tval))))
dimnames(matcoef) <- list(names(tval), c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
cat("\nCoefficient(s):\n")
printCoefmat(matcoef, digits = 6, signif.stars = TRUE)

    cat("\n Degrees of Freedom for the fit:", object$df.fit, "Residual Deg. of Freedom  ", 
        object$df.residual, "\n")
    cat("Global Deviance:    ", format(signif(object$G.deviance)), 
        "\n            AIC:    ", format(signif(object$aic)), "\n            SBC:    ", 
        format(signif(object$sbc)), "\n")
    invisible(object)
}
#####################################################################################
#####################################################################################
